SUBROUTINE DG_RADENE
!-----------------------------------------------------------------------
!
! Definition von DG:   Strahlungsenergiegleichung (0.tes Moment der Strahlungstransportgleichung)
!

      use physco,   only : z1, z2, z3, z4, z12, clight, sqrt2pi, pi, StBolz, grav
      use primvar,  only : X, XZ, XA, DG, MR, MD, ME, MJ, MH, MUr, &
                           DG_i, DG_im1, DG_ip1, DG_im2, DG_ip2
      use global,   only : tst, zz
      use config,   only : np, BCflag
      use geomvar,  only : S_flux, S_vol, S_volZ
      use zvar,     only : rho0, H_p, dH_pdE, dH_pdR
      use advecvar, only : J_adv, dJ_adv_dim2, dJ_adv_dim1, dJ_adv_di, dJ_adv_dip1
      use eddvar,   only : f_edd
      use matvar,   only : Pgas0, dPgas0dE, dPgas0dD, &
                           Tgas, dTgasdE, dTgasdD, &
                           OPAplk, dOPAplkdE, dOPAplkdD


      implicit none

      integer :: i


!-----------------------------------------------------------------------
!    Randbedingungen
!-----------------------------------------------------------------------

   ! innere Pseudozellen: i=1 & i=2
         DG(MJ,MJ,DG_i,1)      =  z1
         DG(MJ,MR,DG_i,1)      =  (X(MD,1)*StBolz*z4 *Tgas(1)**3*dH_pdR(1)*dTgasdD(1)) / (pi*sqrt2pi*H_p(1)**2)
         DG(MJ,MD,DG_i,1)      = -((StBolz*z4 *Tgas(1)**3*dTgasdD(1)) / (pi*sqrt2pi*H_p(1)))
         DG(MJ,ME,DG_i,1)      = -((StBolz*z4 *Tgas(1)**3*(-((X(MD,1)*dTgasdD(1)*dH_pdE(1)) / (sqrt2pi*H_p(1)**2)) + dTgasdE(1))) / pi)

         DG(MJ,MJ,DG_i,2)      =  z1
         DG(MJ,MR,DG_i,2)      =  (X(MD,2)*StBolz*z4 *Tgas(2)**3*dH_pdR(2)*dTgasdD(2)) / (pi*sqrt2pi*H_p(2)**2)
         DG(MJ,MD,DG_i,2)      = -((StBolz*z4 *Tgas(2)**3*dTgasdD(2)) / (pi*sqrt2pi*H_p(2)))
         DG(MJ,ME,DG_i,2)      = -((StBolz*z4 *Tgas(2)**3*(-((X(MD,2)*dTgasdD(2)*dH_pdE(2)) / (sqrt2pi*H_p(2)**2)) + dTgasdE(2))) / pi)

   ! aeussere Pseudozellen + ueberzaehlige skalare Zelle: i=np, i=np-1, i=np-2
         DG(MJ,MJ,DG_i,np)     =  z1
         DG(MJ,MR,DG_i,np)     =  (X(MD,np)*StBolz*z4 *Tgas(np)**3*dH_pdR(np)*dTgasdD(np)) / (pi*sqrt2pi*H_p(np)**2)
         DG(MJ,MD,DG_i,np)     = -((StBolz*z4 *Tgas(np)**3*dTgasdD(np)) / (pi*sqrt2pi*H_p(np)))
         DG(MJ,ME,DG_i,np)     = -((StBolz*z4 *Tgas(np)**3*(-((X(MD,np)*dTgasdD(np)*dH_pdE(np)) / (sqrt2pi*H_p(np)**2)) + dTgasdE(np))) / pi)

         DG(MJ,MJ,DG_i,np-1)   =  z1
         DG(MJ,MR,DG_i,np-1)   =  (X(MD,np-1)*StBolz*z4 *Tgas(np-1)**3*dH_pdR(np-1)*dTgasdD(np-1)) / (pi*sqrt2pi*H_p(np-1)**2)
         DG(MJ,MD,DG_i,np-1)   = -((StBolz*z4 *Tgas(np-1)**3*dTgasdD(np-1)) / (pi*sqrt2pi*H_p(np-1)))
         DG(MJ,ME,DG_i,np-1)   = -((StBolz*z4 *Tgas(np-1)**3*(-((X(MD,np-1)*dTgasdD(np-1)*dH_pdE(np-1)) / (sqrt2pi*H_p(np-1)**2)) + dTgasdE(np-1))) / pi)

         DG(MJ,MJ,DG_i,np-2)   =  z1
         DG(MJ,MR,DG_i,np-2)   =  (X(MD,np-2)*StBolz*z4 *Tgas(np-2)**3*dH_pdR(np-2)*dTgasdD(np-2)) / (pi*sqrt2pi*H_p(np-2)**2)
         DG(MJ,MD,DG_i,np-2)   = -((StBolz*z4 *Tgas(np-2)**3*dTgasdD(np-2)) / (pi*sqrt2pi*H_p(np-2)))
         DG(MJ,ME,DG_i,np-2)   = -((StBolz*z4 *Tgas(np-2)**3*(-((X(MD,np-2)*dTgasdD(np-2)*dH_pdE(np-2)) / (sqrt2pi*H_p(np-2)**2)) + dTgasdE(np-2))) / pi)


!-----------------------------------------------------------------------
!    Restlicher Bereich
!-----------------------------------------------------------------------

      do i=3,np-3

DG(MJ,MR,DG_i,i)=-(X(MJ,i)*pi*(X(MR,i) + X(MR,i+1))) + X(MJ,i)*pi*(X(MR,i+1) - X(MR,i)) - &
clight*X(MH,i)*pi*tst*z2  - f_edd(i)*X(MJ,i)*pi*tst*X(MUr,i)*z2  + (X(MJ,i)*pi*(X(MR,i) + X(MR,i+1))*tst*(X(MR,i)*X(MUr,i) + &
X(MR,i+1)*X(MUr,i+1))*z12 *(-z1 **2 + f_edd(i)*z3 )) / (X(MR,i)**2 + X(MR,i+1)**2) - (X(MJ,i)*S_vol(i)*tst*X(MUr,i)*z12 *(-z1 **2 &
+ f_edd(i)*z3 )) / (X(MR,i)**2 + X(MR,i+1)**2) - (X(MJ,i)*pi*tst*(X(MR,i)*X(MUr,i) + X(MR,i+1)*X(MUr,i+1))*(X(MR,i+1) - &
X(MR,i))*z12 *(-z1 **2 + f_edd(i)*z3 )) / (X(MR,i)**2 + X(MR,i+1)**2) + (X(MJ,i)*X(MR,i)*S_vol(i)*tst*(X(MR,i)*X(MUr,i) + &
X(MR,i+1)*X(MUr,i+1))*z12 *z2 *(-z1 **2 + f_edd(i)*z3 )) / (X(MR,i)**2 + X(MR,i+1)**2)**2 - (-(pi*(XA(MR,i) + X(MR,i))) - &
pi*(X(MR,i) - XA(MR,i)) + pi*tst*XZ(MUr,i)*z2 *zz )*J_adv(i) - clight*pi*rho0(i)*(X(MR,i) + X(MR,i+1))*tst*OPAplk(i)*(X(MJ,i) - &
(StBolz*Tgas(i)**4) / pi) + clight*pi*rho0(i)*tst*(X(MR,i+1) - X(MR,i))*OPAplk(i)*(X(MJ,i) - (StBolz*Tgas(i)**4) / pi) - &
(clight*X(MD,i)*S_vol(i)*tst*OPAplk(i)*(X(MJ,i) - (StBolz*Tgas(i)**4) / pi)*dH_pdR(i)) / (sqrt2pi*H_p(i)**2) - &
(clight*X(MD,i)**2*S_vol(i)*tst*(X(MJ,i) - (StBolz*Tgas(i)**4) / pi)*dH_pdR(i)*dOPAplkdD(i)) / (sqrt2pi**2*H_p(i)**3) + &
(clight*X(MD,i)**2*(X(MR,i) + X(MR,i+1))*StBolz*tst*(X(MR,i+1) - X(MR,i))*z4 *OPAplk(i)*Tgas(i)**3*dH_pdR(i)*dTgasdD(i)) / &
(sqrt2pi**2*H_p(i)**3)
DG(MJ,MR,DG_ip1,i)=X(MJ,i)*pi*(X(MR,i) + X(MR,i+1)) + X(MJ,i)*pi*(X(MR,i+1) - X(MR,i)) + &
clight*X(MH,i+1)*pi*tst*z2  + f_edd(i)*X(MJ,i)*pi*tst*X(MUr,i+1)*z2  - (X(MJ,i)*S_vol(i)*tst*X(MUr,i+1)*z12 *(-z1 **2 + &
f_edd(i)*z3 )) / (X(MR,i)**2 + X(MR,i+1)**2) - (X(MJ,i)*pi*(X(MR,i) + X(MR,i+1))*tst*(X(MR,i)*X(MUr,i) + &
X(MR,i+1)*X(MUr,i+1))*z12 *(-z1 **2 + f_edd(i)*z3 )) / (X(MR,i)**2 + X(MR,i+1)**2) - (X(MJ,i)*pi*tst*(X(MR,i)*X(MUr,i) + &
X(MR,i+1)*X(MUr,i+1))*(X(MR,i+1) - X(MR,i))*z12 *(-z1 **2 + f_edd(i)*z3 )) / (X(MR,i)**2 + X(MR,i+1)**2) + &
(X(MJ,i)*X(MR,i+1)*S_vol(i)*tst*(X(MR,i)*X(MUr,i) + X(MR,i+1)*X(MUr,i+1))*z12 *z2 *(-z1 **2 + f_edd(i)*z3 )) / (X(MR,i)**2 + &
X(MR,i+1)**2)**2 + (-(pi*(XA(MR,i+1) + X(MR,i+1))) - pi*(X(MR,i+1) - XA(MR,i+1)) + pi*tst*XZ(MUr,i+1)*z2 *zz )*J_adv(i+1) + &
clight*pi*rho0(i)*(X(MR,i) + X(MR,i+1))*tst*OPAplk(i)*(X(MJ,i) - (StBolz*Tgas(i)**4) / pi) + clight*pi*rho0(i)*tst*(X(MR,i+1) - &
X(MR,i))*OPAplk(i)*(X(MJ,i) - (StBolz*Tgas(i)**4) / pi)
DG(MJ,MD,DG_i,i)=(clight*S_vol(i)*tst*OPAplk(i)*(X(MJ,i) - (StBolz*Tgas(i)**4) / pi)) / (sqrt2pi*H_p(i)) + &
(clight*X(MD,i)*S_vol(i)*tst*(X(MJ,i) - (StBolz*Tgas(i)**4) / pi)*dOPAplkdD(i)) / (sqrt2pi**2*H_p(i)**2) - &
(clight*X(MD,i)*(X(MR,i) + X(MR,i+1))*StBolz*tst*(X(MR,i+1) - X(MR,i))*z4 *OPAplk(i)*Tgas(i)**3*dTgasdD(i)) / &
(sqrt2pi**2*H_p(i)**2)
DG(MJ,ME,DG_i,i)=-((clight*X(MD,i)*S_vol(i)*tst*OPAplk(i)*(X(MJ,i) - (StBolz*Tgas(i)**4) / pi)*dH_pdE(i)) / &
(sqrt2pi*H_p(i)**2)) + clight*rho0(i)*S_vol(i)*tst*(X(MJ,i) - (StBolz*Tgas(i)**4) / pi)*(-((X(MD,i)*dOPAplkdD(i)*dH_pdE(i)) / &
(sqrt2pi*H_p(i)**2)) + dOPAplkdE(i)) - clight*rho0(i)*(X(MR,i) + X(MR,i+1))*StBolz*tst*(X(MR,i+1) - X(MR,i))*z4 &
*OPAplk(i)*Tgas(i)**3*(-((X(MD,i)*dTgasdD(i)*dH_pdE(i)) / (sqrt2pi*H_p(i)**2)) + dTgasdE(i))
DG(MJ,MUr,DG_i,i)=-(f_edd(i)*X(MJ,i)*pi*X(MR,i)*tst*z2 ) - (X(MJ,i)*X(MR,i)*S_vol(i)*tst*z12 *(-z1 **2 + &
f_edd(i)*z3 )) / (X(MR,i)**2 + X(MR,i+1)**2) - pi*XZ(MR,i)*tst*z2 *zz *J_adv(i)
DG(MJ,MUr,DG_ip1,i)=f_edd(i)*X(MJ,i)*pi*X(MR,i+1)*tst*z2  - (X(MJ,i)*X(MR,i+1)*S_vol(i)*tst*z12 *(-z1 **2 + &
f_edd(i)*z3 )) / (X(MR,i)**2 + X(MR,i+1)**2) + pi*XZ(MR,i+1)*tst*z2 *zz *J_adv(i+1)
DG(MJ,MJ,DG_im2,i)=-(S_flux(i)*zz *dJ_adv_dim2(i))
DG(MJ,MJ,DG_im1,i)=-(S_flux(i)*zz *dJ_adv_dim1(i)) + S_flux(i+1)*zz *dJ_adv_dim2(i+1)
DG(MJ,MJ,DG_i,i)=S_vol(i) + f_edd(i)*pi*tst*(X(MR,i+1)*X(MUr,i+1) - X(MR,i)*X(MUr,i))*z2  - &
(S_vol(i)*tst*(X(MR,i)*X(MUr,i) + X(MR,i+1)*X(MUr,i+1))*z12 *(-z1 **2 + f_edd(i)*z3 )) / (X(MR,i)**2 + X(MR,i+1)**2) + &
clight*rho0(i)*S_vol(i)*tst*OPAplk(i) - S_flux(i)*zz *dJ_adv_di(i) + S_flux(i+1)*zz *dJ_adv_dim1(i+1)
DG(MJ,MJ,DG_ip1,i)=-(S_flux(i)*zz *dJ_adv_dip1(i)) + S_flux(i+1)*zz *dJ_adv_di(i+1)
DG(MJ,MJ,DG_ip2,i)=S_flux(i+1)*zz *dJ_adv_dip1(i+1)
DG(MJ,MH,DG_i,i)=-(clight*pi*X(MR,i)*tst*z2 )
DG(MJ,MH,DG_ip1,i)=clight*pi*X(MR,i+1)*tst*z2

      end do


END SUBROUTINE DG_RADENE
